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We calculate in this article the CMB and matter power spectra of a class of early f(R) cosmologies, 
which takes the form of f(R) = R + \iHoexp[R/(\2Ho )]. Unlike the late-time f(R) cosmologies 
such as f(R) = R + a(-R) 13 (/3 < 0), the deviation from ACDM of this model occurs at a higher 
redshift (thus the name Early f(R) Cosmology) , and this important feature leads to rather different 
ISW effect and CMB spectrum. The matter power spectrum of this model is, at the same time, again 
very sensitive to the chosen parameters, and LSS observations such as SDSS should constrain the 
parameter space stringently. We expect that our results are applicable at least qualitatively to other 
models that produce f(R) modification to GR at earlier times (e.g., redshifts O(10) < z < 0(1)) 
than when dark energy begins to dominate - such models are strongly disfavored by data on CMB 
and matter power spectra. 
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I. INTRODUCTION 

The theory of f(R) gravity has attracted a lot of at- 
tention recently. Theoretically, there are no principles 
prohibiting the inclusion of higher order curvature in- 
variants, such as R 2 , R^R^v and R a ^ l5 Ra/3-yS and their 
scalar functions, into the Einstcin-Hilbert action, as long 
as the various constraints from local and cosmological ob- 
servations are satisfied. Practically, appropriate choices 
of these correction terms could modify the late behav- 
ior of the Universe drastically, thus giving the recently- 
observed accelerating expansion an explanation from an 
alternative theory of gravity. In 0, U the authors con- 
sider a specific model in which the correction is a polyno- 
mial of R 2 , R^R^ and R a ^ s R a p~ l& , and their analysis 
shows that there exist late-time accelerating attractor so- 
lutions. Meanwhile, models with R 2 , R^R^ corrections 
are discussed within the Palatini approach 0, 0, IE El- 
in which the field equations are second order, and simi- 
lar acceleration solutions are found. This Palatini-/(_R) 
theory of gravitation is then tested using various cosmo- 
logical data such as Supernovae (SN), Cosmic Microwave 
Background (CMB) shift parameter, baryon oscillation 
and Big Bang Nucleosynthesis (BBN) , in |2L |g, III - 

Recently, the constraints from CMB and matter power 
spectra are obtained for the typical case of f(R) = 
R + a(-R) 13 in 0,[n|, which successfully exclude most 
of the parameter space, making the model indistinguish- 
able from the standard ACDM in practice. As the cor- 
rections to General Relativity (GR) in this model, which 
act as an effective dark energy component, become im- 
portant only at very late times (or very low energy densi- 
ties) , these works almost exclude the possibility of modi- 
fying theory of gravity in the /(i?)-manner at late times. 
However, there remains the possibility that f(R) mod- 



ification to GR occurs in the earlier time. Obviously, 
the successful predictions of BBN and primary CMB re- 
quire that effects from this modification cannot be sig- 
nificant during these earlier times. Thus we will concen- 
trate in this work on the case where the f(R) modifica- 
tion becomes important during some intermediate times 
between recombination and present. To be specific, we 
shall consider a class of f{R) gravity model with the 
form f(R) = R + \iH§exp(R/ \ 2 Hq), the so-called ex- 
ponential f(R) gravity model, which, as shown below in 
Sec. Ill A, well possesses the desirable property of mod- 
ifying gravitational theory at intermediate cosmic times 
(e nerg y densities). A similar model has been considered 
in |12fl using the metric approach, and an early-time In- 
tegrated Sachs- Wolfe (ISW) effect-Large Scale Structure 
(LSS) cross correlation is found there, which is poten- 
tially useful for distinguishing this model from others. 
Here, however, we shall analyze it using the Palatini ap- 
proach and exactly calculate its CMB and matter power 
spectra by numerically solving the covariant perturba- 
tion equations given in [Tll | ; comparisons with the model 
f(R) = R + a(-R)& and ACDM would also be presented 
and discussed. This work therefore contributes to filling 
in the gap in constraining f(R) cosmology throughout 
the cosmic history. 

This work is organized as the following: in Sec. II we 
briefly introduce the model and list the background and 
perturbation evolution equations that are needed in the 
numerical calculation. In Sec. Ill we incorporate these 
equations into the public code CAMB [13| and give the 
numerical results. The discussion and conclusion are then 
presented in Sec. IV. Throughout this work we will adopt 
the unit c = 1 and only consider the case of a spatially 
flat (background) Universe without massive neutrinos. 



II. FIELD EQUATIONS IN f(R) GRAVITY 
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In this section we briefly summarize the main ingredi- 
ents of f(R) gravity in Palatini formulation and list the 
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general perturbation equations in this theory. 



A. The Modified Einstein Equations 

The starting point of the Palatini-/(i?) gravity is the 
Einstein-Hilbert action 



S = 



d A x, 



-9 



in 



(1) 



in which k = 8ttG, R = g ab R ab {T) (a, b = 0, 1, 2, 3) with 
^ob(r) being defined as 



Rab 



ab,c 



1 ac,b ' 1 cd L ab 1 ad L cb ■ 



(2) 



holds with respect to it, we shall treat this metric as the 
physical one. Consequently the difference between f(R) 
gravity and GR could be understood as a change of the 
manner in which the spacetime curvature and thus the 
physical Ricci tensor lZ ab responds to the distribution 
of matter (through the modified Einstein equations) . In 
order to make this point explicit, we rewrite Eq. (3) with 
the aid of Eq. (7) as 

K%b = FTZab - -jQab! 

+ ^V a FV b F - V a V b F - X -g ah V c V c F, (8) 
in which F, f are now F(T)J{T). 



Note that the connection f is not the Levi-Civita connec- 
tion of the metric g ab , which we denote as T; rather it is 
treated as an independent field in the Palatini approach. 
Correspondingly, the tensor R ab and scalar R are not the 
ones calculated from g ab as in GR, which are denoted by 
lZ a b and 1Z respectively in this work (1Z = g ab lZ ab )- The 
matter Lagrangian density C m , however, is assumed to 
be independent of F, which is the same as in GR. 

Varying the action with respect to the metric g ab then 
leads to the modified Einstein equations 
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FR a b — 7;9abf(R) — K-T ab , 



(3) 



where F = F(R) = df(R)/dR and T ab is the energy- 
momentum tensor. The trace of Eq. (3) reads 



FR — 2/ = kT 



(4) 



with T = p — 3p. This is the so-called structural equation 
which relates R directly to the energy components in 
the Universe: given the specific form of f(R) and thus 
F(R), R as a function of T could be obtained by numer- 
ically or analytically solving this equation. 

The extremization of the action Eq. (1) with respect 
to the connection field f gives another equation 



V a [F(R)J—gg bc ] 



0, 



(5) 



indicating that the connection T is compatible with the 
metric "f a b which is conformal to g ab : 



j ab = F(R)g ab . 



(6) 



With Eq. (6) we could easily obtain the relation between 
R ab and TZ ab 

_ , 3V a FV b F V a V b F g ab V c V c F 

i^ab — 'Vob H 7T^9 ~ A') 



2F 2 



F 



2F 



Note that in the above we use T> and V to denote the 
covariant derivative operators compatible with g ab and 
7 afc respectively. 

Since C m depends only on g ab (and some matter fields) 
and the conservation law of the energy-momentum tensor 



B. The Perturbation Equations 

The perturbation equations in general theories of f(R) 
gravity have been derived in 1 141 . However, here we adopt 
a different derivation by u s llll| which uses the method 
of 1 + 3 decomposition |13l fl5| . For more details we refer 
the reader to [ll| and here we shall simply list the results. 

The main idea of 1 + 3 decomposition is to make space- 
time splits of physical quantities with respect to the 4- 
velocity u a of an observer. The projection tensor h ab is 
defined as h ab = g ab — u a u b which can be used to obtain 
covariant tensors perpendicular to u. For example, the 
covariant spatial derivative T> of a tensor field T^.'^ is 
defined as 



b--c — 
d--e — 



■ ■ h%h r d ■ ■ kv'V:* 



(9) 



The energy-momentum tensor and covariant derivative 
of velocity could be decomposed respectively as 

%b = TTafc + %<l{ a U b ) + pU a U b - ph abl (10) 

V a u b = a ab + w ab + -0h ab + u a A b . (11) 

In the above Tt ab is the projected symmetric trace free 
(PSTF) anisotropic stress, q the vector heat flux, p the 
isotropic pressure, o ab the PSTF shear tensor, vo ab — 
^[aUb], 9 = T^ Cu c = 3d/ a (a is the cosmic scale factor) 
the expansion scalar and A b = ii b the acceleration; the 
overdot denotes time derivative expressed as <f> — u a T> a (j), 
and the square brackets mean antisymmetrization and 
parentheses symmetrization. The normalization is cho- 
sen as u 2 = 1. 

Decomposing the Riemann tensor and making use of 
the modified Einstein equations, we obtain, after lin- 
earization, five constraint equations 



F 







Kq a 



V c (e ab cd u d m ab ); 
3FV a F , 0V a F 



2F 2 

\v a e 



3F 



VgF 

F 



V a ab + V w ab 



(12) 



(13) 
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B ab 



1 



e 6)ec 



2F 



V b ir. 



ab 



2 . F \ 2-a 



1 



-—n{p + p)V a F- 



V b B ab = —k 



D c q d + {p + p)w cd 



cd n ,b 



(14) 

(15) 
(16) 



Here e abcd is the covariant permutation tensor, £ ab and 
B ab are respectively the electric and magnetic parts of 
the Weyl tensor. 

Meanwhile, we obtain seven propagation equations: 

p+{p + p)9 + V a q a = 0;(17) 
q a + ^8q a + (p + p)A a -V aP + f> b ir ab = 0;(18) 



+ 



3F 
2F 



3F 2 3_F np J_ 

2F 2 ~ 2F ~ ~F ~ 2F 



- V a A a 
V 2 F 



2F 



= 0;(19) 



3F 
AF 



Sab + ^ K7T ab + ^pP{aPb) F 

2 

ru+ -8zu - T> [a A b ] 



2F 



0;(20) 
0;(21) 







3P\ 


K a b + 


6> 


-T^yab 



2F 



(p + p)a ab +V {a q b) 







£ah + 


H) 





^d{a e b)ec ue 



B ab +\e + —\B ab + V c £ d{a e b)ec d u e 



0;(22) 



1 

2F 



KV c TT d(a e b)ec d u e 



0,(23) 



where the angle bracket means taking the trace free part 
of a quantity. 

Besides the above equations, it would also be useful to 
express the projected Ricci scalar 1Z in the hypersurfaces 
orthogonal to u a as .11] 



n 



«(P + 3p) - / 



3F 
2F 



2V 2 F 
F 



(24) 



The spatial derivative of the projected Ricci scalar, rj a 
ha'D a 'R-, is then given as 



2F 



F 



—F 

2F 



V a F 



a ^. „ cl ~ , » o _,. 2a 



a 

~3F 







'±F-9 


—F + 9 






2F 




2F 



2F 
V a F, 



F + 9 



v a e 

(25) 



and its propagation equation 
26» 



Va + Y Va 



a 
2F 

a 



F 3 



V a V 2 F~ -nV a V c q c 



F 29 
F + Y 



aV a V c A c .(26) 



Since we are considering a spatially flat Universe, 
the spatial curvature must vanish for large scales which 
means that 7Z = 0. Thus from Eq. (24) we obtain 



F 
2F 



1 2 



= g^[K(p + 3p)-/]. (27) 



This is just the modified (first) Friedmann equation of 
f(R) gravity, and the other modified background equa- 
tions (the second Friedmann equation and the energy- 
conservation equation) could be obtained by taking the 
zero-order parts of Eqs. (17, 19). It is easy to check that 
when f(R) = R, we have F — 1 and these equations 
reduce to those in GR. 

Recall that we have had /, F and R as functions of 
T] it is then straightforward to calculate F, F, T> a F, T> a F 
etc. as functions of T = —(p b + p c )9 and T> a T — 
(1 — 3c 2 )T> a p b + T> a p c , in which p b ^ is the energy density 
of baryons (cold dark matter) and c s is the baryon sound 
speed. Notice that we choose to neglect the small baryon 
pressure except in the terms where its spatial derivative 
is involved, as they might be significant at small scales. 
The above equations could then be numerically propa- 
gated given the initial conditions, telling us information 
about the evolutions of small density perturbations and 
the CMB anisotropics in theories of f(R) gravity. These 
results will be given in the next section. 



III. NUMERICAL RESULTS 

Now let us consider the model f(R) = R + 
AiiJoCxp(i?/A2-ffo )■ There are two positive dimension- 
less parameters in this model, Ai and A2, of which, 
roughly speaking, Ai controls the overall significance of 
the correction to GR while A2 determines the time at 
which the correction becomes important: at early times 
when \R\ ^> \2H 2 , the correction \iHQexp(R/ \ 2 Hq) — > 
and at late times when \R\ <C it tends to a con- 

stant XiHq, thus retrieving the ACDM model; deviation 



from ACDM becomes significant when \R\ ~ ^H^. Solar 
system constraint on A2 is estimated in [I3 to be A 2 <C 
10 6 , which is also applicable here. Ho is a constant with 
dimension Mpc -1 which helps to make Ai, A2 dimension- 
less and we will take it to be the currently preferred value 
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of the Hubble constant, i.e., Hq = 72 km s 1 Mpc 1 , for 
later convenience. 

Just like in the case of f{R) = R + a{— R) 13 , it could 
be easily checked that the parameters Ai, A2 and f2 m (the 
dimensionless density parameter of baryon + cold dark 
matter) are not independent. We thus can choose A2, f2 m 
as the independent parameters and numerically calculate 
Ai from Eqs. (4, 27). For purpose of illustration we will 
fix Q m = 0.3 in this work. 



A. The Background Evolution 

The background evolution of this model with different 
values of A2 is shown in Fig. 1. The results are rather 
similar to those obtained with the metric approach [l2j . 
As expected, the deviation from ACDM is maximized at 
some intermediate ages and reduced towards both earlier 
and later times. At z < 2, the deviations are less than 
- 0.2% and ~ 0.5% for A 2 = 1000 and 500; such small de- 
viations are difficult to be identified by the measurements 
of H(z) such as SNe observations. However, as shown be- 
low, their effects on the CMB and matter power spectra 
could be rather large and A2 = 1000 could be safely ex- 
cluded by current observation on matter power spectrum 
alone. 

The evolution of F(R) — 1 provides another measure- 
ment of the deviation from ACDM because in the latter 
F is simply 1. It is clear that F increases at some ear- 
lier times and finally tends to another constant larger 
than 1. Since, from Eq. (27), the effective gravitational 
constant is G/F, the model becomes effectively ACDM 
with a smaller Newton's constant at its later evolution- 
ary stage. However, we expect this effect to be small 0] 
for large enough A 2 . 



B. The TT CMB Spectrum 

We have calculated the CMB spectrum of the present 
model using CAMB, with adiabatic initial conditions be- 
cause the corrections to GR are expected not to influence 
initial conditions. The results are shown in Fig. 2, where 
we have plotted the CMB spectra for A 2 = 00, 1000 and 
500. Interestingly we see that the deviation from ACDM 
(A2 = 00) occurs roughly in the range I e (10,500), 
while for even smaller Z's it vanishes. This is in con- 
trast to the f{R) = R + a(—R)P model, in which dras- 
tic changes happen at I < 300 |ll|. The contributions 
from ISW effect alone are also shown in Fig. 2, both 
for the f(R) = R + X 1 H^exp(R/ \ 2 H§) and for the 
f(R) = R + a(-R)' 3 models. Obviously these agree qual- 
itatively with the behaviors of the modified CMB spectra 
in the two models, since at low Z's a significant late-time 
contribution comes from the ISW effect. 

The ISW effect involves the integral of the time- 
variation of the potential weighted by the spherical Bessel 




FIG. 1: (Color online) The background evolutions. Upper 
panel: the ratio between the expansion rates in this exponen- 
tial f(R) cosmology and in ACDM as a function of the scale 
factor a. Lower panel: F(R) — 1 as a function of a. The black, 
red, green, blue, cyan and magenta curves represent cases of 
A 2 = 00, 1000, 500, 300, 100, 50 respectively. 

functions for each I 

I{ sw (k) = 2 / e-^' k ji[k(vo-v)]dV, (28) 
Jo 

in which r is the optical depth, rj is the conformal time, a 
prime denotes derivative with respect to rj and subscript 
'0' means the current value; <f>k is Weyl tensor variable 
defined as the coefficient of harmonic expanding £ a b in 
terms of Q k b [H 

/ — 1 a 1 

k 

where Q k b = §2-£>(aA>)Q fc and Q k is the eigenfunction of 
T> 2 with eigenvalue k 2 /a 2 . 

The evolutions of the potential <f>k for various /c's are 
plotted in Fig. 3, and these provide a qualitative expla- 
nation for Fig. 2. For the ISW effect, we know that the 
photons typically travel through many peaks and val- 
leys of the perturbation, and this leads to cancellations 
of the effects. Mathematically this arises from the os- 
cillatory property of the spherical Bessel function ji(x): 
fdxj,(x) = V^r[±(Z + l)]/2r[£(Z + 2)] a (tt/20 1 / 2 M- 
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FIG. 2: (Color online) Upper panel: The theoretical CMB 
spectrum of the present model. From bottom to top: 
A2 = 00, 1000, 500. Lower panel: Contributions from the 
ISW effect. Solid curves are the results of the f(R) = 
R + Ai_ffgexp(ii/A2-ffo ) model, from bottom to top: A2 = 
00, 1000, 500; Dashed curves are those of the f{R) = R + 
a(-R) 13 model for comparison, from bottom to top: j3 = 
0, —0.05, —0.1. In both cases we use Q m — 0.3. 
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FIG. 3: (Color online) Upper panel: The potential (j>k 
as a function of the scale factor for the different fc-modes 
as indicated besides the curves for the f(R) ~ R + 
\iHgexp(R/ \2II0) model. The black, red and green curves 
represent cases of A2 = 00, 1000 and 500 respectively. Lower 
panel: same as above, but for the f(R) = R + a(— R) 13 model, 
and the black, red and green curves here represent cases of 
j3 = 0, —0.01 and -0.05. In both cases we use Q m = 0.3. 



Thus if <p' k is a constant, then I^ sw (k) would roughly 
scale as l~ x / 2 for a specified k. If <j>' k is not a constant, 
on the other hand, it must be weighted by ji (x) at differ- 
ent times in the integration, and, as ji (x) strongly peaks 
at x ~ I, it is weighted more at x ~ I for an /. To be 
more explicit, for the scales we are interested in the ISW 
integral takes the following form 0] 



e Vfcji[fc(?7o - rj)]dr) 



r 4>' k (Vk 



Vi r[i(f + i)] 
1 2k r[i(/ + 2)] : 



in which r) k = t]q — (I + 1 /2) /k is the position of the main 
peak of ji . 

From Fig. 3 it is clear that 4>k decays rapidly earlier 
in the f(R) = R + \iHQexp(R/ \ 2 H§) model than in the 
f(R) = R + a{-Rf model (for such a comparison we 
could choose A2 and (3 such that their endpoints of the 
evolution paths of (f>k are nearly the same in Fig. 3, e.g., 
X-2 = 500 and f3 — —0.01 there). On the other hand, 
at very late times the evolution of <$>\. in the exponential 
f(R) gravity is nearly the same as in ACDM. Now con- 
sider the contribution from a specified fc-mode: if I is very 



large then rjk will be so small that 0' fc (??fe) — 0; thus the 
contribution to this I is negligible. If I is very small, then 
the contribution to this I will be similar to that in ACDM 
as (ji'kiVk) — < / ) feACDM( 7 ?fe)- F° r t ne intermediate l's, since 
the larger I is the smaller 77^. is, some sufficiently larger Vs 
would receive significant contribution from the ISW effect 
if the massive decay of the potential occurs earlier (at a 
time closer to rjk)- Qualitatively this is why in Fig. 2 the 
ISW effect for the f(R) = R+\ 1 H$exp{R/\ 2 H§) model 
extends to larger l's than that of the f(R) = R + a(—R)P 
model and the differences from ACDM are reduced dra- 
matically towards very-small Vs in the case of exponential 
f(R) gravity. 



As for the CMB polarization spectra of the model, 
since they are not as sensitive to the parameters as the 
temperature spectrum, we shall not discuss them here. 
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C. The Matter Power Spectrum 
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FIG. 4: (Color online) The matter power spectrum in the 
present model. The black, red, green, blue and cyan curves 



FIG. 5: (Color online) The matter power spectra at different 
redshifts. From bottom to top: z — 0, z — 5, z — 10 and 
z — 15. The black, red, and green curves represent the cases 
of f(R) =R+ AiH 2 exp(7?/A 2 Ho) (A 2 = 1000), f(R) = R + 
a(-R) 13 (/3 = -0.01) and ACDM respectively. 



As discussed extensively in 0, 0] , the response of this 
modified gravity to the spatial variations of matter dis- 
tribution could be very sensitive (this could also be seen 
from Eq. (24) above), and at small enough scales this 
feature would significantly affect the growth of density 
perturbations. To be explicit, at large enough fc's, the 
growth equation for the como ving energy density fluctu- 
ations 5 m could be written as [lfj 

d 2 S k 2 F' 

"" • - (30) 



dx 2 



H 2 3F(2FH 



F>) 

in which x = log a and F' /3F(2FH + F') — c 2 acts as an 
effective sound speed squared. For the present exponen- 
tial f(R) theory we have seen from Fig. 1 that F' > 0, so 
1 that c 2 > and the small scale density fluctuations will 
oscillate instead of growing according to Eq. (30). On the 
other hand, if c 2 < 0, then these small scale fluctuations 
will become unstable and blow up dramatically |18j ; this 
does happen in the f{R) = R + a(-R) 13 model if [3 is 
positive |ll| , but not for our choices of parameters in the 
f(R) = R + Ai£/o 2 exp(i?/ \ 2 H 2 ) model. 

We have also calculated the matter power spectrum 
in the present model for different choices of A2 with the 
modified CAMB code and the results are shown in Fig. 4. 
As expected, the matter spectrum is more sensitive to the 
deviation from GR than the CMB spectrum, because it 
relates directly to the density perturbations. From this 
figure we see that A2 5 x 10 4 will hopefully be ex- 
cluded when tested against the measurements on matter 
power spectrum from Sloan Digital Sky Survey (SDSS) 
|19| , which is more stringent than the possible CMB con- 
straints. 

It will also be interesting to consider when the devia- 
tions of the matter power spectrum from the ACDM case 
occur in the f(R) gravity theories. Since c 2 = in the 
strict ACDM model, the deviation will appear when F' 
is significantly nonzero, and it thus appears earlier if F 
begins to evolve earlier. Fig. 5 confirms this expectation: 
it is apparent that for the f(R) — R+XiHQexp(R/ X2H 2 ) 
model the deviation from the ACDM matter power spec- 
trum has almost shaped up before z ~ 5 while for the 
f(R) = R + a(-R) 13 model this happens much later. 



IV. DISCUSSION AND CONCLUSIONS 

To summarize, we have investigated the background 
evolution, CMB and matter power spectra of a class 
of early f(R) gravity model, which takes the form of 
f(R) = R+\iHqCxp(R/\ 2 Hq). Such a model has the in- 
teresting feature that the deviation from ACDM happens 
only at some intermediate time, and this is why its CMB 
and matter power spectra are rather different from those 
in either the ACDM case or some late-time f(R) gravity 
theories, such as f(R) = R+a(-Rf. For the CMB spec- 
trum, the earlier massive decay of the potential extends 
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the ISW effect to larger i's, leaving the ISW at very low Vs 
indistinguishable from that in ACDM. Consequently the 
deviation from the latter is most significant at intermedi- 
ate I's (I £ (10,400 - 500)), as shown in Fig. 2. For the 
matter power spectrum, the modifications to GR lead to 
the appearance of effective pressure fluctuations, which 
restricts the growths of small-scale density perturbations 
and leads to oscillations in small scales of the spectrum. 
It is interesting that, the deviation from ACDM matter 
power spectrum develops rather early in the early f(R) 
gravity model because F here begins to evolve earlier. 

Since the matter power spectrum directly shows the 
density perturbations while the CMB spectrum does not, 
the former is more sensitive to corrections to GR in /(i?) 
gravity, as shown in ■ The constraints from local (so- 
lar system) measurements plus the SDSS data might be 
able to restrict the parameter A2 into a rather narrow 
range, roughly 5 ~ 10 x 10 4 < A2 <C 10 6 , but an exact 
constraint obtainable through carefully searching the pa- 



rameter space with, say, the Markov Chain Monte Carlo 
method is beyond the scope of the work here. 

The present work therefore contributes to filling in 
the gap in constraining f(R) corrections to GR at some 
higher redshifts than when the late-time accelerating ex- 
pansion of our Universe begins; it shows that such early 
f(R) cosmologies are strongly disfavored by current data 
on CMB and matter power spectra, just as their late-time 
correspondences such as the model f(R) = R + a(-R) 13 . 
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